Analyse données qPCR - Modélisation Supervisée
B. Govaerts & M. Martin
mars 04, 2018,22:08
1 Description et objectif du cas d’étude
La transcriptomique est l’étude des ARNm qui sont transcrits à partir d’un génome. Elle permet de quantifier le niveau d’expression d’un gène transcrit en ARNm qui est ensuite traduit en protéine. La mesure de l’expression des gènes peut se faire par la méthode de la PCR quantitative (quantitative Polymerase Chain Reaction ou qPCR). On obtient ainsi une mesure relative de l’expression d’un gène par rapport à un contrôle interne (normalisation par rapport à un/plusieurs gène(s) de référence).
Les données à analyser correspondent à des niveaux d’expression de micro-ARN (miR), soit des courtes chaines d’ARN (en moyenne 22 nucléotides) qui sont des régulateurs traductionnels: en s’appariant à un ARNm complémentaire, il va le dégrader ou réprimer sa traduction en protéine.
Des données d’expression de 19 miRs mesurées par qPCR sont disponibles et l’on vous demande de les analyser à l’aide de certaines méthodes multivariées vues durant les sessions théoriques.
Vous disposez pour cela du jeu de données data_qPCR.csv qui comporte les variables “Groupe” et “Identifiant” caractérisant chaque observation, ainsi que les données d’expression de chaque miR.
2 Instructions
- Analysez les données log-transformées avec les méthodes suivantes:
- 1.1. Régression logistique
- 1.2. Régression logistique avec sélection de variables stepwise-forward
- 1.3. Régression logistique avec Ridge
- 1.4. Régression logistique avec Lasso
- 1.5. PLS-DA
- 1.6. OPLS-DA
- 1.7. Random Forest
Avant de faire ces analyses, créez des matrices vides pour stocker efficacement les différentes sorties suivantes pour chaque classifieur (le classifieur est mis en ligne dans la matrice): * Une matrice reprenant le RMSE, RMSE-CV et la dimensionalité du modèle * Une matrice des coefficients de chaque méthode (sauf pour le Random Forest) * Une matrice des valeurs prédites par cross-validation pour chaque observation
Conseil: créez une fonction de cross-validation de type k-fold pour calculer plus facilement le RMSE-CV.
- Comparez des qualités prédictives des différents classifieurs :
- 2.1. Imprimez le tableau reprenant les RMSE(-CV) et la dimensionalité des modèles.
- 2.2. Réalisez un graphique des coefficients obtenus pour chaque classifieur et comparez les graphiques entre eux.
- 2.3. Représentez graphiquement les valeurs prédites par cross-validation à l’aide de line plots légendés et d’un scatterplot.
- 2.4. Dessinez des courbes ROC présentes sur un même graphe avec légende.
- 2.5. Calculez le cut-off optimal de chaque classifieur.
- 2.6. Enfin, présentez une table de confusion pour chaque classifieur avec le cut-off optimal et calculez le pourcentage de précision (accuracy) de ces classifieurs.
3 Importation et data management
# chemin d'accès
data.path <- "/Users/bgovaerts/Dropbox/Partage_OMICS_COURSE/Data_stat2340/Data_qPCR"
# lecture du data frame
dataqPCR <- read.csv(file.path(data.path,"data_qPCR.csv"),sep=';',header=T)
cat("dimensions \n")## dimensionsdim(dataqPCR)## [1] 120 40cat("colnames \n")## colnamescolnames(dataqPCR)## [1] "Groupe" "Identifiant" "miR_001_5p" "miR_002_5p"
## [5] "miR_003_3p" "miR_004_5p" "miR_005_5p" "miR_006_5p"
## [9] "miR_006_3p" "miR_007_5p" "miR_008_5p" "miR_009_5p"
## [13] "miR_010_3p" "miR_011_3p" "miR_012_5p" "miR_013_5p"
## [17] "let_01_5p" "let_02_5p" "miR_014_5p" "miR_015"
## [21] "miR_016" "Log_miR_001_5p" "log_miR_002_5p" "log_miR_003_3p"
## [25] "log_miR_004_5p" "log_miR_005_5p" "log_miR_006_5p" "log_miR_006_3p"
## [29] "log_miR_007_5p" "log_miR_008_5p" "log_miR_009_5p" "log_miR_010_3p"
## [33] "log_miR_011_3p" "log_miR_012_5p" "log_miR_013_5p" "log_let_01_5p"
## [37] "log_let_02_5p" "log_miR_014_5p" "log_miR_015" "log_miR_016"# matice d'expression
X <- as.matrix(dataqPCR[,3:21]) # matrice d'expression
rownames(X) <- dataqPCR[,2]
obsnames=rownames(X)
X <- X[,sort(colnames(X))]
miR_names <- colnames(X)
Xcat=T # indique que les variables sont des catégories bien spécifiques et pas des longueurs d'onde par exemple. Cela a son importance dans les graphiques de coefficients et loadings
Xlog <- as.matrix(dataqPCR[,22:40])
rownames(Xlog) <- dataqPCR[,2]
Xlog <- Xlog[,sort(colnames(Xlog))]
n = dim(X)[1]
m = dim(X)[2]
classY <- rep(0, n)
classY[dataqPCR[,1]=="malade"] = 1
cat("class factor levels \n")## class factor levelsunique(classY)## [1] 1 0table(classY, dataqPCR[,1])##
## classY malade temoin
## 0 0 57
## 1 63 0# modification de matrices
# centrage miRs
XC <- scale(Xlog, center = T, scale = F)
XC <- data.frame(XC)
# On centre Y pour le PLS et l'O-PLS
YC <- scale(classY, center = T, scale = F)4 Estimation des modèles de classification
Nous allons ajuster 6 types de modèles (regression logistique SW, PLS, O-PLS et logisitique ridge et logistique LASSO, random forest) afin de prédire la classe en fonction des MIR.
Pour chacun des modèles, nous allons stocker le RMSE classique (training dataset = test dataset) ainsi que le RMSE obtenu par cross-validation k-fold. Pour certains modèles, nous allons également stocker les valeurs des coefficients.
RMSE <- function(y_true, y_fit) {
n <- length(y_true)
sqrt((1/n)*sum((y_true-y_fit)^2))
}
##################
# kfoldCV
##################
# fonction qui fait de la k-fold CV sur une fonction de régression prédéfinie prédisant y (FUN)
kfoldCV=function(X, Y, k, FUN = ypred_FUN, seed=100, ...) {
# Inputs
# X, la matrice des régresseurs
# Y, le vecteur des réponses
# k, le nombre de segments et
# FUN : la fonction de régression qui prédit Y
# Outputs :
# ypred.cv, les prédictions par CV
# 1. definition de k segments aleatoires
n=dim(X)[1]
set.seed <- seed
indx <- sample(1:n,n)
nseg <- k
sn <- ceiling(n/k) # taille des segments
# 2. boucle sur les sous-echantillons
# pour calculer y-fit par cv
yfit.cv <- numeric()
for (i in 1:nseg) {
if(i < nseg) {indseg <- indx[(i-1)*sn + (1 : sn)]}
else {indseg <- indx[((i-1) * sn + 1) : n]}
Y.train <- Y[-indseg]
X.train <- X[-indseg,]
X.valid <- X[indseg,]
yfitcv <- FUN(X = X.train, Y = Y.train, XV = X.valid, ...)
yfit.cv[indseg] <- yfitcv
}
# 3. Sortie de la fonction
return(yfit.cv)
}# Définition de paramètres de modélisation
# nombre de segments en cross validation
k <- 10
# Nombre de composantes max
NCmax=10
# Liste des modèles
modnames <- c('Logistique-SW','PLSR','O-PLS','Logistique-RIDGE','Logistique-LASSO','Random Forest')
nmodels <- length(modnames) # nombre de modèles évalués
# Preparation des matrices de résultats
# Matrice des RMSE
RMSEMat <- matrix(nrow = nmodels, ncol = 3)
colnames(RMSEMat) <- c('RMSE-train','RMSE-CV','Taille du modèle')
rownames(RMSEMat) <- modnames
# Matrice des coefficients (pour tous les modèles sauf RF)
CoefficientMat <- matrix(0, ncol = m, nrow = nmodels-1)
colnames(CoefficientMat) <- colnames(XC)[1:m]
rownames(CoefficientMat) <- modnames[-nmodels]
# Matrice des prédictions
PredMat <- matrix(nrow = n, ncol = nmodels)
colnames(PredMat) <- modnames
# Tout est mis dans une liste qui sera remis à jour modèle par modèle
listr=list(CoefficientMat=CoefficientMat,RMSEMat=RMSEMat,PredMat=PredMat)# création d'une fonction pour stocker et imprimer les résultats de chaque modèle
printresmod=function(i,YC,yfit,yfitCV,coef,parmod,listr,pcoef=T)
{
# i = N° du modèle
# YC : réponses observée
# yfit : réponses prédites
# yfitCV : réponses prédites par CV
# coef : vecteur des coefficients
# parmod paramètre de dimension ou de lissage du modèle
# pcoef : T/F pour indiquer si des coefficient sont dispo pour le modèle (pas le cas pour les RF)
par(mfrow=c(1,2))
plot(YC,yfit,pch=20,main="Y-Yfit")
plot(YC,yfitCV,pch=20,main="Y-YfitCV",ylab="yfitCV")
# calcul des RMSE
n=length(YC)
listr$RMSEMat[i,1]=sqrt(sum((YC-yfit)^2)/n)
listr$RMSEMat[i,2] =sqrt(sum((YC-yfitCV)^2)/n)
listr$RMSEMat[i,3] = parmod
# Récupération des coefficients et des réponses
if(pcoef==T) listr$CoefficientMat[i,]=coef
listr$PredMat[,i]=yfit
# affichage de résultats
pander(listr$RMSEMat[i,])
par(mfrow=c(1,1))
if(pcoef==T){
if(Xcat==F){plot(Xval,listr$CoefficientMat[i,],type="l",ylab="Coefficients",main=paste("Coefficients ",dimnames(coefficient)[[1]][i]),xlab=Xxlab)
abline(h=0)}
if(Xcat==T){dotchart(coef,labels=colnames(coef),
xlab="Coefficients",main=paste("Coefficients - ",modnames[i])) }
}
return(listr)
}4.1 (1) Régression logistique stepwise forward
# Ajustement du modèle de regression stepwise forward
null <- glm(classY ~ 1, data=XC,family="binomial") # modele de depart
full <- glm(classY ~ ., data=XC,family="binomial") # modele avec toutes les variables
fit.sw <- step(null, scope=list(lower=null, upper=full), direction="both",trace=0) # critère d'opimisation= AIC
pander(summary.glm(fit.sw))| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.5942 | 0.3373 | -1.761 | 0.07817 |
| log_miR_002_5p | -12.62 | 3.313 | -3.81 | 0.0001387 |
| log_miR_009_5p | 13.01 | 3.304 | 3.939 | 8.197e-05 |
| log_miR_011_3p | 2.022 | 0.9884 | 2.046 | 0.04076 |
| log_miR_010_3p | -1.175 | 0.6565 | -1.79 | 0.07349 |
(Dispersion parameter for binomial family taken to be 1 )
| Null deviance: | 166.1 on 119 degrees of freedom |
| Residual deviance: | 104.1 on 115 degrees of freedom |
yfit.train <- fit.sw$fitted.values
namessw=names(coefficients(fit.sw))
coef <- CoefficientMat[1,]
coef[namessw[-1]]=coefficients(fit.sw)[-1]
nc <- length(coef)-1
ypred_FUN <- function(X,Y,XV,...){
null <- glm(Y ~ 1, data = X,family = "binomial")
full <- glm(Y ~ ., data = X,family = "binomial")
fit.sw <- step(null, scope = list(lower = null,
upper = full),
direction="both",trace=0)
coef <- coefficients(fit.sw)
XV <- cbind(1, XV[,names(coef)[-1]])
yfit.cv <- as.matrix(XV) %*% coef
yfit.cv <- exp(yfit.cv)/(1+exp(yfit.cv))
return(yfit.cv)
}
yfit.cv <- kfoldCV(X = XC, Y = classY, k = k, FUN = ypred_FUN)
# Sauvegardes et impression de résultats
listr=printresmod(1,classY,yfit.train,yfit.cv,coef,nc,listr)
4.2 (2) Partial least squares
PLSR par la librairie pls. On estime directement les modèles par cross validation afin de trouver le bon nombre de composantes. Ensuite on reestime ce modèle là.
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
ypred_FUN <- function(X, Y, XV, ...){
fit.plsr <- plsr(Y ~ . , data = X, ...)
coef <- coefficients(fit.plsr)
yfit.train <- as.matrix(XV) %*% coef
return(yfit.train)
}
## Boucle sur le nombre de composantes possibles: 1 à NCmax
RMSE.cv <- numeric()
yfit.cv <- list()
for(i in 1:NCmax) {
res.cv <- kfoldCV(X = XC, Y = YC, k = k,
FUN = ypred_FUN,
ncomp = i)
yfit.cv[[i]] <- res.cv
RMSE.cv[i] <- RMSE(y_true = classY, y_fit = yfit.cv[[i]])
}
# Recherche du nombre de composantes à garder
plot(1:10, RMSE.cv, type = "o", main = "RMSE cross-validation")
ncomp.opt = which.min(RMSE.cv)
ncomp.opt = max(2, ncomp.opt)
abline(v = ncomp.opt, col = 2, lty = 2)
RMSE.cv <- RMSE.cv[ncomp.opt]
yfit.cv <- yfit.cv[[ncomp.opt]] + mean(classY)
# Ajustement du modèle avec le nombre optimal de composantes
fit.plsr <- plsr(YC ~ . , data = data.frame(XC),
ncomp = ncomp.opt)
print(summary(fit.plsr))## Data: X dimension: 120 19
## Y dimension: 120 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
## 1 comps 2 comps
## X 36.76 47.46
## YC 26.76 33.89
## NULLyfit.train <- fit.plsr$fitted.values[, , ncomp.opt] + mean(classY)
coef <- coefficients(fit.plsr)[, , ]
nc <- ncomp.opt
# Sauvegardes et impression de résultats
listr=printresmod(2,classY,yfit.train,yfit.cv,coef,nc,listr)
4.2.1 Scores
Représentation graphique des scores pour les 2 premières composantes avec la librairie MBX
par(mfrow=c(1,1))
# graphe des scores avec la librairie MBXUCL
ScatterPlot(fit.plsr$scores[,1],fit.plsr$scores[,2], createWindow=FALSE, points_labs =obsnames, main = paste("PLS score plot for PC1 and PC2"), color=classY, pch=classY,xlab="PC1",ylab="PC2")
4.2.2 Graphe des loadings
loadings <- fit.plsr$loadings
for (i in 1:2) {
plot(loadings[,i], type="h", xaxt="n", xlab="",
ylab = "Loading", main = paste0("Loading ", i), lwd = 3)
axis(side = 1, at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)
}
4.3 (3) Orthogonal partial least squares
fit.opls <- opls(x = Xlog, y = classY, predI = 1,
orthoI = 4, scaleC = "center",
printL = TRUE, plotL = FALSE)
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
ypred_FUN=function(X,Y,XV,...) {
fit.opls <- opls(x = X, y = Y, predI = 1,
scaleC = "center" , ...)
ropls::predict(fit.opls, XV)
}
## Boucle sur le nombre de composantes possibles
RMSE.cv <- numeric()
yfit.cv <- list()
for(i in 1:NCmax) {
rescv <- kfoldCV(Xlog, classY, k = k, FUN = ypred_FUN,
orthoI = i, printL = FALSE,
plotL = FALSE)
yfit.cv[[i]] <- rescv
RMSE.cv[i] <- RMSE(y_true = classY, y_fit = yfit.cv[[i]])
}
# Recherche du nombre de composantes à garder
plot(RMSE.cv, type = "o", main = "RMSE cross-validation")
ncomp.opt <- which.min(RMSE.cv)
ncomp.opt <- max(2, ncomp.opt)
abline(v = ncomp.opt, col = 2, lty = 2)
yfit.cv <- yfit.cv[[ncomp.opt]]
# Ajustement du modèle avec le nombre optimal de composantes
fit.opls <- opls(x = Xlog, y = classY, predI = 1, orthoI = ncomp.opt,
scaleC = "center", printL = TRUE, plotL = FALSE)
yfit.train <- ropls::predict(fit.opls, X)
coef <- fit.opls@weightMN %*% fit.opls@cMN
nc <- ncomp.opt
# Sauvegardes et impression de résultats
listr=printresmod(3,classY,yfit.train,yfit.cv,coef,nc,listr)
4.3.1 Scores
Représentation graphique des scores OPLS
scores_p <- fit.opls@scoreMN
scores_o <- fit.opls@orthoScoreMN
par(mfrow=c(1,1))
# graphe des scores avec la librairie MBXUCL
ScatterPlot(scores_p[,1], scores_o[,1], createWindow=FALSE, points_labs =obsnames, main = paste("O-PLS score plot for PCP and PCO1"), color=classY, pch=classY,xlab="PCPredictive",ylab="PCO1")
4.3.2 Graphe des loadings
loadings_p <- fit.opls@loadingMN
loadings_o <- fit.opls@orthoLoadingMN
plot(loadings_p, type="h", xaxt="n", xlab="", ylab = "Loading",
main ="Loading prédictif", lwd = 3)
axis(side = 1, at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)
plot(loadings_o[,1], type="h", xaxt="n", xlab="", ylab = "Loading",
main ="Loading ortho 1", lwd = 3)
axis(side = 1, at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)
4.4 (4) Régression logistique avec RIDGE
# Ajustement du modèle de regression pénalisé avec toutes les données
# lambda2 => RIDGE
fit.ridge <- penalized(response = classY,
penalized = XC,
model = "logistic",
lambda1 = 0,
lambda2 = 0.1,
trace = FALSE)
# Optimisation du paramètre lambda2 par cross-validation
lambda2 <- optL2(response = classY, penalized = XC,
lambda1 = 0, model = "logistic", trace=0)$lambda
# Regression logistique ridge avec lambda optimisé
fit.ridge <- penalized(response = classY,penalized = XC, model = "logistic",
lambda1 = 0, lambda2 = lambda2, trace = FALSE)
yfit.train = attr(fit.ridge,'fitted')
coef <- attr(fit.ridge,'penalized')
nc = sum(!(coef)==0)
# Calcul du yfit.cv
# fonction pour obtenir les y prédits
ypred_FUN = function(X, Y, XV, ...){
fit.penalized <- penalized(response = Y, penalized = X,
model = "logistic", trace=FALSE, ...)
penalizedcoeff <- attr(fit.penalized,'penalized')
intercept <- attr(fit.penalized,'unpenalized')
coef <- c(intercept, penalizedcoeff)
yfit.cv=as.matrix(cbind(1,XV))%*%coef
yfit.cv = exp(yfit.cv)/(1+exp(yfit.cv))
return(yfit.cv)
}
yfit.cv <- kfoldCV(X, classY, k = k, FUN = ypred_FUN, lambda1 = 0,
lambda2 = lambda2)
# Sauvegardes et impression de résultats
listr=printresmod(4,classY,yfit.train,yfit.cv,coef,nc,listr)
4.5 (5) Régression logistique avec LASSO
# Ajustement du modèle de regression pénalisé avec toutes les données
# lambda1=> LASSO
fit.lasso <- penalized(response = classY, penalized = XC, model = "logistic",
lambda1 = 0.1, lambda2 = 0, trace = FALSE)
# Optimisation du paramètre lambda2 par cross-validation
lambda1 <- optL1(response = classY, penalized = XC,
lambda2 = 0, model = "logistic",
trace=0)$lambda
# Regression logistique ridge avec lambda optimisé
fit.lasso <- penalized(response = classY, penalized = XC, model = "logistic",
lambda2 = 0, lambda1 = lambda1, trace = FALSE)
yfit.train = attr(fit.lasso,'fitted')
coef <- attr(fit.lasso,'penalized')
nc = sum(!(coef)==0)
# Calcul du yfit.cv
# fonction pour obtenir les y prédits (c'est la m^ême que la précédente)
ypred_FUN = function(X, Y, XV, ...){
fit.penalized <- penalized(response = Y, penalized = X,
model = "logistic", trace=FALSE, ...)
penalizedcoeff <- attr(fit.penalized,'penalized')
intercept <- attr(fit.penalized,'unpenalized')
coef <- c(intercept, penalizedcoeff)
yfit.cv=as.matrix(cbind(1,XV))%*%coef
yfit.cv = exp(yfit.cv)/(1+exp(yfit.cv))
return(yfit.cv)
}
yfit.cv <- kfoldCV(X, classY, k = k, FUN = ypred_FUN,
lambda1 = lambda1, lambda2 = 0)
# Sauvegardes et impression de résultats
listr=printresmod(5,classY,yfit.train,yfit.cv,coef,nc,listr)
4.6 (6) Random Forest
fit.rf <- randomForest(formula=as.factor(classY) ~ ., data=Xlog, ntree=500, mtry=3)
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
ypred_FUN <- function(X,Y,XV,...) {
fit.rf <- randomForest(formula=as.factor(Y) ~ ., data = X, ntree = 500, ...)
yfit.cv <- predict(fit.rf, newdata = XV, type = 'class')
yfit.cv
}
## Boucle sur le nombre de composantes possibles
RMSE.cv <- numeric()
yfit.cv <- list()
j <- 1
mtry_tested <- seq(3,19, 2)
for(i in mtry_tested) {
rescv <- kfoldCV(Xlog, classY, k = k, FUN = ypred_FUN, mtry = i)
yfit.cv[[j]] <- rescv
yfit <- as.numeric(yfit.cv[[j]]) - 1
RMSE.cv[j] <- RMSE(y_true = classY, y_fit = yfit)
j <- j + 1
}
# Recherche du nombre de composantes à garder
plot(mtry_tested,RMSE.cv, type = "o", main = "RMSE cross-validation")
mtry.opt <- mtry_tested[which.min(RMSE.cv)]
abline(v = mtry.opt, col = 2, lty = 2)
yfit.cv <- as.numeric(yfit.cv[[which.min(RMSE.cv)]])-1
# Ajustement du modèle avec le nombre optimal de mtry
fit.rf <- randomForest(formula=as.factor(classY)
~ ., data = Xlog,
ntree=500, mtry=mtry.opt)
yfit.train <- predict(fit.rf, newdata = Xlog, type='class')
yfit.train <- as.numeric(yfit.train)-1
nc = m
# Sauvegardes et impression de résultats
listr=printresmod(6,classY,yfit.train,yfit.cv,coef,nc,listr,pcoef=F)
5 Evaluation des qualités des classifieurs
5.1 Matrice des RMSE en estimation et k-fold-validation
cat(c("Taille des segments", k))## Taille des segments 10pander(listr$RMSEMat)| RMSE-train | RMSE-CV | Taille du modèle | |
|---|---|---|---|
| Logistique-SW | 0.3829 | 0.443 | 18 |
| PLSR | 0.406 | 0.4301 | 2 |
| O-PLS | 0.5736 | 0.4473 | 2 |
| Logistique-RIDGE | 0.399 | 0.4442 | 19 |
| Logistique-LASSO | 0.4174 | 0.458 | 5 |
| Random Forest | 0 | 0.5477 | 19 |
5.2 Plots des coefficients
for(i in 1:(nmodels-1)) {
dotchart(listr$CoefficientMat[i,],labels=colnames(listr$CoefficientMat),
xlab="Coefficients",main=paste("Coefficients - ",modnames[i]))
} 
5.3 Lines plots of predictions
# Line plot to the Predictions
par(xpd = TRUE, mar = c(2,2,2,7))
plot(1:n,listr$PredMat[,1],lty=1,type="l",
main="Predictions ordered by observation ID",xlab="Observation ID",ylab="Prediction",
ylim=c(min(listr$PredMat),max(listr$PredMat)), col=1)
for(i in 2:nmodels){
lines(1:n,listr$PredMat[,i],lty=i,col=i)
}
legend("topright",1,modnames,col=1:nmodels,lty=1:nmodels,cex=0.8,
inset =c(-0.18,0))
abline(v=92)
Sur ce premier graphique, l’axe des x représente l’ordre dans lequel les mesures ont été effectuées. la droite verticale représente un changement dans le kit utilisé pour mesurer les expresisons des miRs. On peut constater que cela a un impact non négligeable sur les prédictions en cross-validation.
par(xpd = TRUE, mar = c(2,2,2,7))
idx <- order(classY)
plot(1:n,listr$PredMat[idx,1],lty=1,type="l",
main="Predictions ordered by class",xlab="Observation ID",ylab="Prediction",
ylim=c(min(listr$PredMat),max(listr$PredMat)), col=1, xaxt="n")
axis(side=1, at = 1:n, labels = idx, las= 2, cex.axis=0.8)
for(i in 2:nmodels){
lines(1:n,listr$PredMat[idx,i],lty=i,col=i)
}
legend("topright",1,modnames,col=1:nmodels,lty=1:nmodels,cex=0.8, inset =c(-0.18,0))
abline(v=57)
L’axe des x représente un vecteur ordonné de la classe des observations. La ligne verticale délimite les témoins des malades.
5.4 Scatter Plot matrix of predictions:
# Scatter Plot matrix of predictions
pairs(listr$PredMat,main="Scatter Plot Matrix of predictions")
5.5 Courbes ROC
# Courbes de ROC
cutoff <- seq(min(listr$PredMat)-0.1,max(listr$PredMat)+0.1,by=c(0.01))
nco <- length(cutoff)
truepositive <- matrix(nrow=nco,ncol=nmodels)
falsepositive <- matrix(nrow=nco,ncol=nmodels)
npos <- sum(classY==1)
nneg <- sum(classY==0)
listpos <- (classY==1)
listneg <- (classY==0)
for(i in 1:nco){
for(j in 1:nmodels) {
# Rate of true positives in the set of positive
truepositive[i,j] <- sum(listr$PredMat[listpos,j]>=cutoff[i])/npos
# Rate of false positives in the set of negatives
falsepositive[i,j] <- sum(listr$PredMat[listneg,j]>=cutoff[i])/nneg
}
}
plot(falsepositive[,1],truepositive[,1],type="l",main="ROC Curves Q-PCR Models",
xlab="False Positive Rate",ylab="True Positive Rate")
for(j in 2:nmodels){
lines(falsepositive[,j],truepositive[,j],type="l",col=j,lty=j)
}
legend(0.6,0.75,dimnames(listr$PredMat)[[2]],col=1:nmodels,lty=1:nmodels)
5.6 Recherche du cut-off optimal pour chaque méthode
# Recherche des seuils optimaux qui minimisent la distance au point (0,1)
cutoffopt=listr$PredMat[1,]
for(i in 1:nmodels){
distances2 <- falsepositive[,i]^2+(1-truepositive[,i])^2
cutoffopt[i] <- cutoff[which(distances2==min(distances2))][1]
}
pander(t(cutoffopt))| Logistique-SW | PLSR | O-PLS | Logistique-RIDGE | Logistique-LASSO |
|---|---|---|---|---|
| 0.5572 | 0.5472 | 0.6272 | 0.5472 | 0.5472 |
| Random Forest |
|---|
| 0.007205 |
# Dessin des cutoff sur un scatter plot
for(i in 1:nmodels){
plot(YC,listr$PredMat[,i],main=paste("Cut-offs - Method:",
dimnames(listr$PredMat)[[2]][i]),ylab="Prediction")
abline(h=cutoffopt[i])
} 
5.7 Matrices de confusion avec les cut-offs optimaux
PredictClass <- matrix(nrow = n, ncol = nmodels)
dimnames(PredictClass) <- dimnames(listr$PredMat)
TrueClass <- rep(1,n)
for(i in 1:nmodels){
PredictClass[,i] <- listr$PredMat[,i]>=cutoffopt[i]
}
# Calcul les nombres de bien classés avec les seuils optimaux
ClassTable <- matrix("",nrow = nmodels * 3, ncol = 2)
dimnames(ClassTable)[[1]] <- rep(c("","Class0","Class1"),nmodels)
dimnames(ClassTable)[[2]] <- c("FALSE","TRUE")
dimnames(ClassTable)[[1]][((1:nmodels)*3)-2] <- modnames
for(i in 1:nmodels){
ClassTable[(i-1)*3+(2:3),] <- table(classY,PredictClass[,i])
}
cat("Classification Tables")## Classification Tablespander(ClassTable)| FALSE | TRUE | |
|---|---|---|
| Logistique-SW | ||
| Class0 | 42 | 15 |
| Class1 | 12 | 51 |
| PLSR | ||
| Class0 | 39 | 18 |
| Class1 | 11 | 52 |
| O-PLS | ||
| Class0 | 38 | 19 |
| Class1 | 10 | 53 |
| Logistique-RIDGE | ||
| Class0 | 43 | 14 |
| Class1 | 13 | 50 |
| Logistique-LASSO | ||
| Class0 | 42 | 15 |
| Class1 | 16 | 47 |
| Random Forest | ||
| Class0 | 57 | 0 |
| Class1 | 0 | 63 |
# accuracy
pc.accuracy <- c()
for(i in 1:nmodels){
pc.accuracy[i] <- sum(diag(table(classY,PredictClass[,i])))/n
}
pc.accuracy <- as.matrix(pc.accuracy)
colnames(pc.accuracy) <- "pc.accuracy"
rownames(pc.accuracy) <- modnames
pander(pc.accuracy)| pc.accuracy | |
|---|---|
| Logistique-SW | 0.775 |
| PLSR | 0.7583 |
| O-PLS | 0.7583 |
| Logistique-RIDGE | 0.775 |
| Logistique-LASSO | 0.7417 |
| Random Forest | 1 |